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Note: This is a working paper which will be expanded/updated frequently. The directory 
deleeuwpdx.net/pubfolders/isotone has a pdf copy of this article, the complete Rmd hie that 
includes all code chunks, and R and FORTRAN hies with the code. Suggestions are welcome 
24/7. 

1 Problem 

In least squares m.onotone regression with a simple order we minimize a loss function of the 
form 

n 

a (x) = J2 W i( X i ~ Vi) 2 (!) 

i= 1 

over X\ < x 2 < • • • < x n . The vector w contains positive weights. The vector y is the target. 
Monotone regression projects the target on the cone defined by the inequality constraints. 

There are many R implementations available in various GRAN packages, starting with 
isoregO in the stats package (if the weights are all equal). The function amalgmO in this 
paper is another such implementation, merely an R wrapper for a double precision modification 
of the FORTRAN code of Cran (1980). We could also have used the C implementation in recent 
versions of the smacof package. Various additions and improvements of the original package 
of De Leeuw and Mair (2009) are discussed in Mair, De Leeuw, and Groenen (2015). 
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If there are ties in the data the required order on x needs to take those into account. There 
are three basic approaches. The first two are described by Kruskal (1964), the third is in De 
Leeuw (1977). In the primary approach. In the primary approach we require ordering only 
between tie blocks, not within tieblocks. The secondary approach requires equality within tie 
blocks and ordering between tie blocks, the tertiary approach only requires ordering of the 
tie block means. 

So if the data used to generate the order are 

2 2313312 1 
then in the primary approach we require 

{x 4 ,x 7 ,x 9 } < {xi,x 2 ,x 8 } < {x 3 ,x 5 ,x 6 }, 
in the secondary approach this becomes 

X 4 = X 7 = Xg < X\ = x 2 = x 8 < £3 = X5 = Xe, 
and in the tertiary approach it is 

X 4 + X 7 + Xg X! + X 2 + X 8 X 3 + X 5 + X 6 

3 ~ 3 “ 3 

2 Implementation 

The code in the appendix implements the three approaches to ties by using a list of indices, 
computed by 

f <- sort(unique (x)) 

g <- lapply(f, function (z) which(x == z)) 

Thus if the data are 


## 

[1] 2 

1 

3 2 1 3 3 1 2 

then the list 

is 

## 

[[1]] 



## 

[1] 2 

5 

8 

## 




## 

[[2]] 



## 

[1] 1 

4 

9 

## 




## 

[[3]] 



## 

[1] 3 

6 

7 


The function isotoneO uses simple list manipulations to get the data in the correct form 
for a call to amalgmO. Study it at your leisure. 
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3 Example 

Again we use the same data to generate the order 
## [1] 2 1 3 2 1 3 3 1 2 
We also choose the target 
## [1] 7 1 5 6 2 9 3 8 4 
The primary approach gives 

## [1] 5.5 1.0 5.5 5.5 2.0 9.0 5.5 5.5 5.5 

with a loss function value of 17.5. For the secondary approach the projection is 

## [1] 5.666667 3.666667 5.666667 5.666667 3.666667 5.666667 5.666667 3.666667 
## [9] 5.666667 

with a loss function value of 52. And for the tertiary approach we find 
## [1] 7 1 5 6 2 9 3 8 4 

with a loss function value of only 0. We show the transformations of the data in figure 1. 
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Figure 1: Three approaches to ties 

Note that in the tertiary aproach the transformation of the data is generally not monotone. 
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4 Appendix: Code 


amalgm <- function (x, w = rep (1, length (x))) { 
dyn.load ("pava.so") 
n <- length (x) 
a <- rep (0, n) 
b <- rep (0, n) 
y <- rep (0, n) 

If <- .Fortran ("AMALGM", n = as.integer (n), x = as.double (x), w = as.double (w), 
return (lf$y) 

} 

isotone <- 
function (x, 

y> 

w = rep (1, length (x)), 
ties = "secondary") { 
f <- sort(unique (x)) 
g <" lapply(f, function (z) 
which(x == z)) 
n <- length (x) 
k <- length (f) 
if (ties == "secondary") { 
w <- sapply (g, length) 
h <- lapply (g, function (x) 
yW) 

m <- sapply (h, sum) / w 
r <- amalgm (m, w) 
s <- rep (0, n) 
for (i in 1: k) 
sCgCEi]]] <- r[i] 

> 

if (ties == "primary") { 

h <- lapply (g, function (x) 

yW) 

m <- rep (0, n) 
for (i in 1: k) { 

ii <- order (h[ [i]]) 
g[[i] ] <- g[[i]] [ii] 
h[[i]] <- h[[i]] [ii] 

} 

m <- unlist (h) 

r <- amalgm (m, w) 

s <- r [order (unlist (g))] 
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> 

if (ties == "tertiary") { 
w <- sapply (g, length) 
h <- lapply (g, function (x) 

yM) 

m <- sapply (h, sum) / w 
r <- amalgm (m, w) 
s <- rep (0, n) 
for (i in 1: k) 

stgCCi]]] <- y [g [ [i] ] ] + (r [i] - m[i]) 

> 

return (s) 
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